#!/usr/bin/env python3
# -*- coding:utf-8 -*-

'''
author: maoxuner
e-mail: maoxuner@126.com
date: 2016年4月20日
filename: stability_analysis.py
version: 3.2.0
details: 本程序采用的算法为折线滑动面法（滑楔法），
         目测适用于采用滑楔法的各种土石心墙坝的坝坡稳定分析。
         可考虑地震荷载的影响。最终求得Kc的最小值。
         考虑了心墙位置的影响，计算滑坡时，心墙不一起滑动。
         必须配套数据文件 data.py 一起使用，否则程序无法正常运行！
info: 看了VB版本的边坡稳定分析的代码，感觉还是比较简单的。
      但是VB版本的代码并没有考虑心墙的位置，而是把心墙等同于坝壳进行运算的。
      最近在学python，便决定改写成对应版本的，顺便练练手。
warning: 您可以任意复制、传播本程序，但是请保留此注释内容。
         如果您要基于此程序进行二次开发，请补充注释修改的内容，
         并把新的程序无偿的分享给大家。谢谢！
'''

# 坐标范围生成器
def get_range(start, end, step = 1):
    tmp = start
    while tmp <= end:
        yield tmp
        tmp += step

# 判断数据格式是否正确
def is_num(tmp):
    num_items=('1', '2', '3', '4', '5', '6', '7', '8', '9', '0', '.')
    if tmp == '':
        return False
    # 非字符数据时转换成字符数据
    if not isinstance(tmp, str):
        tmp = str(tmp)
    n = 0
    for i in tmp:
        # 两个及以上小数点为假
        if i == '.':
            n += 1
            if n > 1:
                return False
        if i not in num_items:
            return False
    return True

# 输入数据
def get_data(name,data):
    tmp = input('请输入{0}的值：'.format(name))
    # 当输入为空时，保留默认值
    if tmp == '':
        return float(data)
    while not is_num(tmp):
        tmp = input('输入的数据格式不正确！\n请重新输入{0}的值：'.format(name))
    return float(tmp)

# 计算分块重量
def cal_weight(datas, locates):
    w1 = 0.5 * (locates['xd'] - locates['yd'] * datas['m_dam']) * (locates['yd'] - locates['ya']) * datas['v_wet']
    w2 = 0.5 * (locates['xd'] - locates['yd'] * datas['m_dam']) * (locates['ye'] - locates['yd']) * datas['v_nature']
    w3 = 0.5 * (locates['ye'] - locates['yd']) * (locates['xb'] - locates['xd']) * datas['v_nature']
    w4 = 0.5 * (locates['xc'] - locates['xb']) * (locates['yc'] - locates['yd']) * datas['v_nature']
    return w3 + w4, w1 + w2

# 解一元二次方程组
def cal_equation(a, b, c):
    from math import sqrt
    if a == 0:
        if b != 0:
            return - c / b
        elif c != 0:
            return None
    else:
        delta = b ** 2 - 4 * a * c
        if delta < 0:
            return None
        else:
            x1 = (- b + sqrt(delta)) / (2 * a)
            x2 = (- b - sqrt(delta)) / (2 * a)
            return x1, x2

# 计算安全系数（正常/地震）
def cal_factors(datas, locates, weights, f1 = 0, f2 =0):
    from math import sin, cos, tan
    # 两块土体底面的抗剪强度都取坝壳砂砾料平均天然休止角（弧度）
    phi1 = phi2 = datas['m_sand']
    a1 = locates['alpha1']
    a2 = locates['alpha2']
    w1 = weights[0]
    w2 = weights[1]
    a = (w1 * sin(a1) + f1 * cos(a1)) * cos(a1 - a2) + w2 * sin(a2) + f2 * cos(a2)
    b = - w1 * cos(a1) * cos(a1 - a2) * tan(phi1) - (w1 * sin(a1) + f1 * cos(a1)) * sin(a1 - a2) * tan(phi2) - w2 * cos(a2) * tan(phi2) + f1 * sin(a1) * cos(a1 - a2) * tan(phi1) + f2 * sin(a2) * tan(phi2)
    c = w1 * cos(a1) * sin(a1 - a2) * tan(phi1) * tan(phi2) - f1 * sin(a1) * sin(a1 - a2) * tan(phi1) * tan(phi2)
    return a, b, c

# 地震荷载的计算
def earthquake(datas, locates, weights, level):
    factors_quake = {'7': (3.0, 0.1), '8': (2.5, 0.2), '9': (2.0, 0.4)}
    # 根据地震等级选取相关系数
    am, ah = factors_quake[str(level)]
    # 地震作用的效应折减系数
    XI = 0.25
    # 两块土体底面的抗剪强度都取坝壳砂砾料平均天然休止角
    w1 = weights[0]
    w2 = weights[1]
    # 求荷载作用点
    def center(y1, y2, y3, y4 = 0):
        if y4 == 0:
            return (y1 + y2 + y3) / 3
        else:
            return (y1 + y2 + y3 + y4) / 4
    # 求作用点的动态系数
    def cal_ai(hi, h):
        if h <= 40:
            return 1 + hi / h * (am -1)
        elif hi <= 0.6 * h:
            return 1 + hi / (0.6 * h) * (am - 1) / 3
        else:
            return 1 + (am - 1) / 3 + (hi - 0.6 * h) / (0.4 * h) * (am - (1 + (am - 1) / 3))
    h2 = center(locates['ya'], locates['yd'], locates['ye'])
    if locates['xb'] == locates['xc']:
        h1 = center(locates['yd'], locates['ye'], locates['yb'])
    else:
        h1 = center(locates['yd'], locates['ye'], locates['yb'], locates['yc'])
    f1 = ah * XI * w1 * cal_ai(h1, datas['h_dam'])
    f2 = ah * XI * w2 * cal_ai(h2, datas['h_dam'])
    return f1, f2

def get_location(datas):
    from math import atan
    from data import STEP
    yb = datas['h_dam']
    yc = yb
    xb = yb * datas['m_dam']
    for h_drop in get_range(datas['h_drop'], datas['h_water'], STEP):
        yd = h_drop
        x_core = xb + 0.5 * datas['w_dam'] - 0.5 * datas['w_core']
        ya_max = yd
        xd_min = yd * datas['m_dam'] + STEP
        if datas['h_water'] < datas['h_core']:
            x_water_core = x_core - (datas['h_core'] - datas['h_water']) * datas['m_core']
            xd_max = min(xb, x_water_core)
        else:
            xd_max = xb
        for ya in get_range(0, ya_max, STEP):
            xa = ya * datas['m_dam']
            for xd in get_range(xd_min, xd_max, STEP):
                alpha2 = atan((yd - ya) / (xd - xa))
                xe = xd
                ye = xd / datas['m_dam']
                xc_min = xb
                xc_max_core = xd + (x_core - xd) / (datas['h_core'] - yd) * (yc - yd)
                if ya == yd:
                    xc_max_trend = xb + datas['w_dam']
                else:
                    xc_max_trend = xd + (xd - xa) / (yd - ya) * (yc - yd)
                xc_max = min(xb + datas['w_dam'], xc_max_core, xc_max_trend)
                for xc in get_range(xc_min, xc_max, STEP):
                    alpha1 = atan((yc - yd) / (xc - xd))
                    yield xa, ya, xb, yb, xc, yc, xd, yd, xe, ye, alpha1, alpha2

# 主程序，定义变量，调用程序，输出结果
def main():
    # 从数据文件 data.py 中导入基础数据，请确保 data.py 跟主程序处于同一目录
    try:
        from data import keys_datas, contents, datas, keys_locates, locates, kc, level
    except ImportError:
        print('\n数据文件 data.py 丢失！')
        print('请从 http://git.oschina.net/maoxuner/stability_analysis 下载最新的版本！')
        exit()
    print('请按照以下提示内容来输入数据，不用手动修改 data.py 中的内容')
    # 判断是否要输入数据
    if input('是否使用预置数据（y/n）') != 'y':
        for i in keys_datas:
            datas[i] = get_data(contents[i], datas[i])
        level = int(get_data('地震级别（0,7~9）', level))
    # 输出程序运行时所使用的数据
    print('请确认您本次计算所使用的数据，如果有错误，请重新运行本程序！')
    for i in range(len(keys_datas)):
        print(contents[keys_datas[i]], '：', datas[keys_datas[i]])
    print('地震级别：', level)
    if level not in (7, 8, 9):
        print('地震等级不在程序计算范围内，本次计算将不考虑地震作用。')
        level = 0
    with open('out.txt', 'w') as f:
        # 生成坐标
        for tmp in get_location(datas):
            for i in range(len(keys_locates)):
                locates[keys_locates[i]] = tmp[i]
                f.write('%12.8f,' % tmp[i])
            weights = cal_weight(datas, locates)
            if level == 0:
                factors = cal_factors(datas, locates, weights)
            else:
                forces = earthquake(datas, locates, weights, level)
                factors = cal_factors(datas, locates, weights, forces[0], forces[1])
            result = cal_equation(factors[0], factors[1], factors[2])
            f.write('%12.8f,%12.8f\n' % result)
            if result[0] > 0 and result[0] < kc:
                kc = result[0]
                final = locates.copy()
        # 输出数据
        for i in range(len(keys_datas)):
            f.write('\n{0}：{1}'.format(contents[keys_datas[i]], datas[keys_datas[i]]))
        f.write('\n地震级别：{0}'.format(level))
        print('\n安全系数为：%10.8f\n对应各点的坐标为:' % kc)
        f.write('\n安全系数为：%10.8f\n对应各点的坐标为:\n' % kc)
        for i in range(len(keys_locates)):
            print('%6s %8.4f ' % (keys_locates[i], final[keys_locates[i]]), end = '')
            f.write('%6s %8.4f ' % (keys_locates[i], final[keys_locates[i]]))
            if i % 2 == 1:
                print()
                f.write('\n')

if __name__ == '__main__':
    print('Hello, World!')
    print('技术为懒人而生！欢迎使用本程序！')
    print('本程序采用的算法为折线滑动面法（滑楔法），目测适用于采用滑楔法的各种土石心墙坝。')
    flag = '6'
    while flag != 'q':
        main()
        flag = input('直接回车将继续计算，输入q则退出程序：')

